Основная сложность работы с такими рядами, как в нашем проекте, заключается в том, как учесть сложную структуру сезонности — суточной, недельной и годовой. В рамках моделей ARIMA можно учесть только одну из них. Обычно в таких случаях сезонность с самым маленьким периодом явно моделируют с помощью аримы, а все остальные учитывают за счёт регрессионной компоненты.
Для учёта недельной сезонности используются регрессионные признаки следующего вида: $$s_i = \sin \left( [1,\dots,T] * 2 \pi i/ 168 \right), c_i = \cos \left( [1,\dots,T] * 2 \pi i/ 168 \right), i=1,\dots, K.$$ Здесь $T$ — это длина моделируемого ряда, 168 — длительность недели в часах, а значение параметра $K$ вам предстоит подобрать самостоятельно (в зависимости от длины подготовленного ряда, можно для начала взять K равным 2-5).
Если вы собрали данные за несколько лет, аналогичные признаки можно использовать для учёта годовой сезонности. Длина года — 8766 часов.
Если в данных есть линейный или описываемый ещё какой-то простой функцией $f$ тренд, стоит добавить к регрессионным признакам вектор $[1,…,T]$ или, соответственно, #f([1,…,T])#.
Если вы забыли, как в statsmodels работать с моделями ARIMA, почитайте туториал от создателей.
Чтобы сдать задание, выполните следующую последовательность действий.
import pandas as pd
import numpy as np
from scipy import stats
import statsmodels.api as sm
import seaborn as sns
import statsmodels.formula.api as smf
import statsmodels.stats.api as sms
import scipy as sc
from sklearn import metrics
import itertools
import os
from glob import glob
import pickle
from tqdm import tqdm_notebook
import warnings
warnings.filterwarnings('ignore')
%pylab inline
%%time
# загружаем данные (за 6 месяцев)
f_lst = glob('agr_tripdata_*.pkl')
data = pd.DataFrame()
for f in f_lst:
data = data.append(pd.read_pickle(f))
print(data.shape)
data.info()
# выбираем для анализа самую "нагруженную" ячейку
cell_idx = data.groupby('bins').sum().idxmax()
print('выбираем ячейку {0}, из которой за период совершено {1:,.0f} поездок'
.format(cell_idx[0], data.groupby('bins').sum().max()[0]))
# выбираем данные по выбранной ячейки
cell_data = pd.DataFrame(data.numb[data.bins == cell_idx[0]].values,
index=data.date[data.bins == cell_idx[0]].values,
columns=['numb'])
cell_data.sort_index(inplace=True)
cell_data.head()
cell_data.info()
# рассмотрим ряд изначальный
cell_data.plot(figsize=(16,8), grid=True, legend=False)
plt.ylabel('amount')
pylab.show()
"""
Ярко выражена сезонность суточная и недельная. Дисперсия, будем считать, в стабилизации не нуждается
Ярко выраженного тренда тоже не наблюдаю. Есть пара провалов: Рождество и аномальный снегопад в Нью-Йорке 23.01.2016
И в конце мая - праздник День памяти (30.05.2016) после двух выходных
Были попытки сделать отдельные признаки по праздникам, но как-то я не нашел эту затею дельной... Может не прав
"""
pass
plt.figure(figsize(18,12))
sm.tsa.seasonal_decompose(cell_data).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(cell_data.numb/60)[1])
"""
Критерий Дики-Фуллера сразу отвергает гипотезу нестационарности. Но все не просто, как и обещали.
Трендовая компонента имеет явную недельную сезонность. Будем пытаться от нее избавиться с помощью регрессии
"""
pass
# процедура формирования массива признаков
def gen_features(cell_data, ii):
# генерируем значения
sin_i = [np.sin(np.array(range(1, len(cell_data)+1))*2*np.pi*i/168)
for i in ii]
cos_i = [np.cos(np.array(range(1, len(cell_data)+1))*2*np.pi*i/168)
for i in ii]
# наименования переменны x_[_s, _c] - синус/косинус [_i] - множитель
name_s = ['x_s_'+str(i) for i in ii]
name_c = ['x_c_'+str(i) for i in ii]
# соберем в таблицу, добавим отклик и проиндексируем
features = pd.DataFrame(sin_i + cos_i).T
features.columns = name_s + name_c
features['y'] = cell_data.numb.values
features.index = cell_data.numb.index
return features
# формируем признаки
ii = range(2,7)
features = gen_features(cell_data, ii)
"""
При включении в список семерки SARIMAX падает по исключению singular matrix. Поэтому ограничился шестеркой
Смотрел и боле высокие частоты, что только удаляло от цели - очистить остатки от недельной сезонности
"""
features.head(3)
# добавим бинарный признак выходных
features['weekend'] = features.index.weekday
features['weekend'] = (features['weekend'] < 4) # сначала выделял выходные, потом инвертировал (не переименовал)
# от признаков Рождества и периода снежного шторма отказался
features['holidays'] = features.index.isin(pd.date_range('2015-12-25', periods=10*24, freq='h'))
features['anomality'] = features.index.isin(pd.date_range('2016-01-23', periods=2*24, freq='h'))
# посмотрим на результат на интервале в 4 недели
fig, axes = plt.subplots(figsize=(18, features.shape[1]*2), nrows = features.shape[1], sharex=True)
for idx, name in enumerate(features.columns):
axes[idx].plot(features[name][:168*4], label=name)
axes[idx].legend()
plt.subplots_adjust(hspace=0)
%%time
# воспользуемся statmodels
m1 = smf.ols('y~'+'+'.join(features.columns.drop(['y', 'x_c_3', 'x_c_5'])),
data=features)
fitted = m1.fit(cov_type='HC1')
print(fitted.summary())
plt.figure(figsize(16,7))
plt.subplot(121)
sc.stats.probplot(fitted.resid, dist="norm", plot=pylab)
plt.subplot(122)
fitted.resid.plot.hist()
plt.xlabel('Residuals')
plt.figure(figsize(16,5))
fitted.resid.plot()
plt.title('residuals')
plt.figure(figsize(16,5))
fitted.fittedvalues.plot()
plt.title('fitted values')
pylab.show()
"""
'x_c_3' и 'x_c_5' алгоритм счел статистически незначимыми
Модель имеет небольшую склонность к занижению прогноза. Причем чаще это максимальное занижение.
Завышение имеет больший разброс, и чаще величина завышения небольшая. В будние дни модель завышает чаще.
Из-за этого визуально в остатках не удалось избавиться от недельной сезонности (
"""
pass
# посмотрим на результаты на двухнедельном интервале
fig, axes = plt.subplots(figsize=(18,5))
plt.plot(features.y[:168*4]);
plt.title('target')
fig, axes = plt.subplots(figsize=(18,5))
plt.plot(fitted.resid[:168*4]);
plt.title('residuals');
fig, axes = plt.subplots(figsize=(18,5))
plt.plot(fitted.fittedvalues[:168*2]);
plt.title('predictions');
"""
Выходные стали чуть менее резко отличаться от будних дней (сравните target и residuals).
Так что чуть сгладить недельную сезонность все-таки получилось
"""
pass
# построим STL-декомпозицию остатков
plt.figure(figsize(18,12))
sm.tsa.seasonal_decompose(fitted.resid).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(fitted.resid)[1]);
"""
Критерий Дики-Фуллера отвергает гипотезу нестационарности ряда.
Но в трендовой компоненте мы видим явную недельную сезонность
Будем дифференцировать!
"""
pass
residuals = pd.DataFrame(fitted.resid, columns=['resid'])
residuals.head(3)
Здесь момент творчества. Есть два варианта сезонных лагов 168 и 24. Использовал в разных комбинациях, остановился на этом
# сначала дифференцируем с сезонным лагом 24
residuals['resid_diff_s'] = residuals.resid - residuals.resid.shift(24)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(residuals.resid_diff_s.dropna()).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(residuals.resid_diff_s.dropna())[1]);
"""
В трендовой компоненте осталася явная недельная сезонность
"""
pass
# теперь дифференцируем с сезонным лагом 168
residuals['resid_diff_s_1'] = residuals.resid_diff_s - residuals.resid_diff_s.shift(168)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(residuals.resid_diff_s_1.dropna()).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(residuals.resid_diff_s_1.dropna())[1])
"""
Трендовой компонента выглядит уже лучше...
"""
pass
# продифференцируем еще с единичным лагом
residuals['resid_diff_s_2'] = residuals.resid_diff_s_1 - residuals.resid_diff_s_1.shift(1)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(residuals.resid_diff_s_2.dropna()).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(residuals.resid_diff_s_2.dropna())[1])
"""
Трендовой компонента все более похожа на "случайную пульсацию". Есть всплески, но надо уже остановиться,
т.к. многократное дифференцирование - не есть хорошо
"""
pass
# подбираем параметры модели
plt.figure(figsize(18,10))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(residuals.resid_diff_s_2.dropna().values.squeeze(), lags=168, ax=ax, zero=False)
pylab.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(residuals.resid_diff_s_2.dropna().values.squeeze(), lags=168, ax=ax, zero=False)
pylab.show()
"""
Дифференцирование остатков, которым страдали выше, было ради выбора параметров.
Но в выборе параметров есть еще одно существенное ограничение - терпение в ожидании окончания расчета
Тот факт, что SARIMAX бывает падает через часок расчетов по исключению singular matrix (не часто),
не добавляет желания увеличивать количество параметров. Поэтому на коррелограммы я, конечно, посмотрел,
но гиперпараметры выбрал по своему усмотрению, чуть ниже рекомендуемых
"""
pass
Начальные приближения установим следующие:
ps = range(0, 9)
d = 1 # одно простое
qs = range(0, 6)
Ps = range(0, 1)
D = 2 # два сезонных
Qs = range(0, 1)
parameters = itertools.product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
%%time
# заготовка из C5W1
results = []
best_aic = float("inf")
# c tqdm не так утомительно ждать; если не установлен, переставьте комментарий
#for param in parameters_list:
for param in tqdm_notebook(parameters_list):
#try except нужен, потому что на некоторых наборах параметров модель не обучается
try:
model=sm.tsa.statespace.SARIMAX(cell_data,
order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], 24),
exog=features[features.columns.drop(
['y', 'x_c_3', 'x_c_5'])].astype({'weekend':'int8',
'holidays':'int8',
'anomality':'int8'})
).fit(disp=-1)
#выводим параметры, на которых модель не обучается и переходим к следующему набору
except ValueError:
print('wrong parameters:', param)
continue
aic = model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True).head())
# разобьем график помесячно, так лучше видно
st = datetime.datetime.strptime("2015-12-01", "%Y-%m-%d")
for i in range(1,7):
plt.figure(figsize(18,5))
cell_data[st:st+relativedelta(months=1)].plot(label='data', legend=True)
best_model.fittedvalues[st:st+relativedelta(months=1)].plot(label='predict', legend=True)
st = st+relativedelta(months=1);
print(best_model.summary())
plt.figure(figsize(15,8))
plt.subplot(211)
best_model.resid.plot()
plt.ylabel(u'Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(best_model.resid.values.squeeze(), lags=48, ax=ax, zero=False)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(best_model.resid, 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(best_model.resid[168+24*2:])[1])
print('среднечасовое значение поездок из выбранной зоны - {0:.2f}, корень из среднеквадратичной ошибки прогноза - {1:.2f}'
.format(cell_data.mean()[0], sqrt(metrics.mean_squared_error(cell_data, best_model.fittedvalues))))
Выводы:
остатки не смещены - подверждается критерием Стьюдента (не отвергнута гипотеза о равенстве среднего нулю);
остатки стационарны - подтверждается критерием Дики-Фуллера (отвергнута гипотеза о нестационарности остатков) и визуально
остатки автокоррелированы - подтверждается критерием Льюнга-Бокса (отвергается гипотеза о том, что данные являются "белым шумом") и коррелограммой - имеются статистически значимые отклонения от нуля). Особенно смущает сильно автокоррелированный сезонный лаг. Что-то есть в данных, что не учитывает наша модель.
Дисперсия остатков великовата. Соотношение среднечасового количества поездок, и среднеквадратичной ошибки прогноза указывает на то, что процентов на 20% мы в среднем ошибаемся. Многовато... Но, с другой стороны, мы не мажем совсем мимо! То есть прогноз не выглядит абсолютным бредом :)
Попытки улучшить прогноз на участках праздников и прочих аномалий (добавление признаков), существенных улучшений не дали.
Возможно это все следствие волюнтаризма в выборе гиперпараметров, но времени для экспериментов оказалось недостаточно